Nonlinear and Perturbative Evolution of Distorted Black Holes. II. Odd-parity Modes 
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We compare the fully nonlinear and perturbative evolution of nonrotating black holes with odd- 
parity distortions utilizing the perturbative results to interpret the nonlinear results. This introduc- 
tion of the second polarization (odd-parity) mode of the system, and the systematic use of combined 
techniques brings us closer to the goal of studying more complicated systems like distorted, rotating 
black holes, such as those formed in the final inspiral stage of two black holes. The nonlinear evolu- 
tions are performed with the 3D parallel code for Numerical Relativity, Cactus, and an independent 
axisymmetric code, Magor. The linearized calculation is performed in two ways: (a) We treat the 
system as a metric perturbation on Schwarzschild, using the Regge- Wheeler equation to obtain the 
waveforms produced, (b) We treat the system as a curvature perturbation of a Kerr black hole 
(but here restricted to the case of vanishing rotation parameter a) and evolve it with the Teukolsky 
equation The comparisons of the waveforms obtained show an excellent agreement in all cases. 
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I. INTRODUCTION 

Coalescing black holes are considered one of the most 
promising sources of gravitational waves for gravitational 
wave observatories like the LIGO/VIRGO/GEO/TAMA 
network under construction (see, e.g., Ref. and ref- 
erences therein). Reliable waveform information about 
the merger of coalescing black holes can be crucial not 
only to the interpretation of such observations, but also 
could greatly enhance the detection rate. Therefore, it 
is crucial to have a detailed theoretical understanding of 
the coalescence process. 

It is generally expected that full scale, 3D numerical 
relativity will be required to provide such detailed infor- 
mation. However, numerical simulations of black holes 
have proved very difficult. Even in axisymmetry, where 
coordinate systems are adapted to the geometry of the 
black holes, black hole systems are difficult to evolve be- 
yond about t = 150M, where M is the mass of the sys- 
tem Q. In 3D, the huge memory requirements, and in- 
stabilities associated presumably with the formulations 
of the equations themselves, make these problems even 
more severe. The most advanced 3D calculations based 
on traditional Cauchy evolution methods published to 
date, utilizing massively parallel computers, have diffi- 
culty evolving Schwarzschild Q, Misner or distorted 
Schwarzschild beyond about t = 50M. Characteris- 
tic evolution methods have been used to evolve distorted 
black holes in 3D indefinitely ||^, although it is not clear 
whether the technique will be able to handle highly dis- 
torted or colliding black holes due to potential trouble 



with caustics. 

In spite of such difficulties, much physics has been 
learned and progress has been made in black hole simu- 
lations, in both axisymmetry and in 3D. In axisymmetry, 
calculations of distorted black holes with |l^ and with- 
out angular momentum [ pl]p2| , Misner two black hole 
initial data, including variations of boosted and unequal 
mass black holes, have been all been successfully carried 
out |p|, p^|Jl4|] , and the waveforms generated during the 
collision process have been extensively compared to cal- 
culations performed using perturbation theory [p"5[-p^. 
In 3D, similar calculations have been carried out, espe- 
cially in evolution of 3D distorted black holes |^ where it 
was shown that very accurate waveforms can be extracted 
as a distorted black holes settles down, as is expected to 
happen when two black holes coalesce. 

One of the important results to emerge from these 
studies is that the full scale numerical and perturbative 
results agree very well in the appropriate regimes, giving 
great confidence in both approaches. In particular, the 
perturbative approach turned out to work extremely well 
in some regimes where it was not, a priori, expected to 
be accurate. For example, in the head-on collision of two 
black holes (using Misner data), the perturbative results 
for both waveforms and energy radiated turned out to be 
remarkably accurate against full numerical simulations - 
even in some cases where the black holes had distinct 
apparent horizons. These impressive agreements have 
then been improved by the use of second order pertur- 
bation theory (see Ref. for a comprehensive review 
on the Zerilli approach and Ref. EOl for the more recent 
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curvature based approach holding also for rotating black 
holes). Study of perturbations also offered the plausi- 
ble explanation that the peak of the potential barrier 
that surrounds a black hole was the more relevant quan- 
tity, not the horizon. In a more complex application, 
the collision of boosted black holes was studied. With a 
small boost, the total energy radiated in the collision was 
shown to go down when compared with the Misner data. 
Linear theory was able to show that there were two com- 
ponents to the radiation, one from the background Mis- 
ner geometry and one from the boost. These two compo- 
nents are anti-correlated and combined to produce what 
has since been called the "Baker dip" . Had there only 
been a perturbative analysis one might have worried that 
nonlinear effects might eliminate the dip. Had there been 
only a full numerical simulation, the dip might have been 
thought to be evidence of a coding error. When the two 
were combined, however, a confirmation of the correct- 
ness of both procedures was established and the effect 
understood. 

These are just two examples of a rather large body of 
work that has led to a revival of perturbative calcula- 
tions, now considered to be used as a tool to aid in the 
verification and interpretation of numerically generated 
results. The potential uses of this synergistic approach 
to black hole hole evolutions, combining both numerical 
and perturbative evolutions, are many. First, the two 
approaches go hand-in-hand to verify the full scale non- 
linear numerical evolutions, which will become more and 
more difhcult as 3D binary mergers of unequal mass black 
holes are attempted, with linear, spin, and orbital angu- 
lar momentum. Second, as the above examples show, 
they can aid greatly in the interpretation and physical 
understanding of the numerical results, as also shown in 
3D distorted black hole evolutions . Such insight will 
become more important as we move towards more com- 
plex simulations. (As an example of this below, we will 
show how nonlinear effects and mode-mixing can be un- 
derstood and cleanly separated from linear effects with 
this approach.) 

Finally, there are at least two important ways in which 
a perturbative treatment can actually aid the numerical 
simulation. First, as shown in Ref. it is possible to 
use perturbative evolutions to provide good outer bound- 
ary conditions for a numerical simulation, since away 
from the strong field region one expects to see low am- 
plitude gravitational waves propagating on a black hole 
background. This information can be exploited in the 
outer region in providing boundary data. Second, this 
combined approach can be used in future applications of 
perturbative approaches to "take over" and continue a 
previously computed full scale nonlinear numerical sim- 
ulation. For example, if gravitational waveforms are of 
primary interest in a simulation, once the system has 
evolved towards a perturbative regime (e.g., two coa- 
lescing black holes form a distorted Kerr hole, or evolve 
close enough that a close limit approximation is valid), 
then one may be able to extract the relevant gravitational 



wave data, and evolve them on the appropriate black hole 
background to extract waveforms |22[| . Not only would 
such a procedure save computational time, it may actu- 
ally be necessary in some cases to extend the simulations. 
As discussed above, 3D black hole evolutions using tra- 
ditional ADM style formulations, with singularity avoid- 
ing slicings, generally break down before complete wave 
forms can be extracted. A perturbative approach may be 
necessary in such cases to extract the relevant waveform 
physics. This work (called Lazarus project) is currently 
being undertaken by some of the present authors. 

However, all work to date in this area of comparing full 
scale numerical simulations with perturbative approaches 
has dealt with even-parity distortions of Schwarzschild- 
like black holes. See for instance Ref. , referred to here 
as Paper I, where we compared perturbative techniques, 
based on the Zerilli approach, with fully nonlinear evolu- 
tions of even-parity distorted black holes. This restriction 
to the Zerilli equation cannot handle the odd-parity class 
of perturbations, and more importantly, it cannot be ap- 
plied easily to the case of rotating black holes. The more 
general black hole case has both even- and odrf-parity 
distortions, and also involves black holes with angular 
momentum. For this reason, in this paper we take an 
important step towards application to the more general 
case of rotating, distorted black holes, by introducing 
the Teukolsky equation as the fundamental perturbation 
equation. In fact, for black holes with angular momen- 
tum, there is not an i — m multipole decomposition of 
metric perturbations in the time domain and the most 
natural way to proceed is with the curvature-based per- 
turbation formalism leading to the Teukolsky equation, 
which also simultaneously handles, in a completely gauge 
invariant way, both even- and odd-parity perturbations. 

The paper is structured as follows. In section II, we 
review the initial data sets and four different techniques 
and approaches to evolve black holes. 

1. We first carry forward the metric-based perturba- 
tion approach by considering the Regge- Wheeler 
(odd-parity) equation to perform perturbative evo- 
lutions, and for the first time apply these tech- 
niques to a class of distorted black hole data sets 
containing even-and odd-parity distortions. 

2. We also show how one can carry out such perturba- 
tive evolutions with the curvature-based Teukolsky 
equation, using the same initial datasets. Although 
in certain cases, the metric perturbations can be 
computed from the curvature perturbations, and 
vice versa psf , in general using both approaches 
helps us better to understand the systems we are 
dealing with. 

3. We carry out fully nonlinear evolutions of the same 
data sets for comparison with a 2D (axisymmet- 
ric) code, Magor, also capable of evolving distorted 
rotating black holes. 
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4. Finally, the same initial data is evolved in its full 
3D mode with a general parallel code for numerical 
relativity, Cactus. 

In section III we finally discuss the results in detail and 
show how the combination of these different approaches 
provides an extremely good and systematic strategy to 
cross check and further verify the accuracy of the codes 
used. The comparisons of waveforms obtained in this 
way show an excellent agreement, in both perturbative 
and full nonlinear regimes. 

Although in this paper we restrict ourselves to the case 
of initial datasets without angular momentum, the family 
of datasets we use for this study also includes distorted 
Kerr black holes, which will be considered in a follow 
up paper. In fact, our eventual goal is to apply both 
fully nonlinear numerical and perturbative techniques to 
evolve a binary black hole system near the merger phase, 
which final stage can be reasonable modeled by a single 
distorted Kerr black hole. In this case we should be able 
to address extremely important questions like how much 
energy and angular momentum can be radiated in the 
final merger stage of two black holes. 



II. FOUR WAYS TO EVOLVE DISTORTED 
BLACK HOLES 

A. Distorted Black Hole Initial Data 

Our starting point is represented by a distorted black 
hole initial data sets developed originally by Brandt and 
Seidel [p|-|lOt to mimic the coalescence process. These 
data sets correspond to "arbitrarily" distorted rotating 
single black holes, such as those that will be formed in the 
coalescence of two black holes. Although this black hole 
family can include rotation, in this first step we restrict 
ourselves to the non-rotating limit (the so-called "Odd- 
Parity Distorted Schwarzschild" of Ref. How- 
ever, these data sets do include both degrees of gravita- 
tional wave freedom, including the "rotation-like" odd- 
parity modes. 

The details of this initial data procedure are covered 
in [p|-|Tot, so we will go over them only briefly here. We 
follow the standard 3+1 ADM decomposition of the Ein- 
stein equations which give us a spatial metric, an extrin- 
sic curvature, a lapse and a shift. We choose our system 
such that we have a conformally flat three-metric 7^ de- 
fined by 



* = ( 1 



(1) 



where the coordinates and are the usual spherical 
coordinates and the radial coordinate has been replaced 
by an exponential radial coordinate?] (f = ^e"^)- Thus, 
if we let the conformal factor be ^ = ^/f we have the flat 
space metric with the origin at 77 = —00. If we let 
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we have the Schwarzschild 3-metric. In this case one finds 
that 77 = corresponds to the throat of a Schwarzschild 
wormhole, rj = ±00 corresponds to spatial infinity in each 
of the two spaces connected by the Einstein-Rosen bridge 
(wormhole) . Note also that this metric is invariant under 
the isometry operation 77 — > —rj. In the full nonlinear 2D 
evolution will use this fact to give ourselves the appro- 
priate boundary conditions for making distorted black 
holes. 

The extrinsic curvature is chosen to be 



He 
He He 



where 

HE=qG{{n' + l)-{'2 + n')s\n^ 
He = —djjqc cos 9 sin" 



QG = Qq 



exp 



+ exp (^-~ (t] + riof /a' 



(3) 

(4) 
(5) 

(6) 



The various functions have been chosen so that the mo- 
mentum constraints are automatically satisfied, and have 
the form of odd-parity distortions in the black hole ex- 
trinsic curvature. The function qQ provides an adjustable 
distortion function, which satisfies the isometry opera- 
tion, and whose amplitude is controlled by the parameter 
Qo- This parameter carries units of length squared. Since 
we will be comparing cases with different masses we will 
refer to an amplitude Qo = Qo/-^^ normalized by the 
ADM mass of the initial slice. If Qo vanishes, an unper- 
turbed Schwarzschild black hole results. The parameter 
n' is used to describe an "odd-parity" distortion. It must 
be odd, and have a value of at least 3. The function is 
the conformal factor, which we have abstracted from the 
metric and extrinsic curvature according to the factor- 
ization given by Lichncrowicz p4| . This decomposition 
is valuable, because it allows us to solve the momentum 
and Hamiltonian constraints separately (with this factor- 
ization the extrinsic curvature given above analytically 
solves the momentum constraints). 

For the class of data considered here the only non- 
trivial component of the momentum constraints is the ip 
component: 



d^Hs sin^ + de {He sin^ 6^ 



0. 



(7) 



Note that this equation is independent of the function qQ. 
This enables us to choose the solutions to these equations 
independently of our choice of metric perturbation. 

At this stage we solve the Hamiltonian constraint nu- 
merically to obtain the appropriate value for 4*. Data at 
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the inner boundary {t] = 0) is provided by an isometry 
condition, namely, that the metric should not be changed 
by an inversion through the throat described by ry ^ —77. 
If we allow Qq to be zero, we recover the Schwarzschild 
solution for 5*. 

The Hamiltonian constraint equation can be expanded 
in coordinate form to yield, in this case. 



Hi sin^ e + H] 



(8) 



This construction is similar to that given in Bowen and 
York |2^, except that form of the extrinsic curvature is 
different. The same procedure described above can also 
be used to construct Kerr and distorted Kerr black holes, 
as described in Ref. ]^-[lc|], but we defer that application 
to a future paper. 

Note that although the form of the extrinsic curva- 
ture is decidedly odd-parity (consider reversing the the 
(^-direction), the Hamiltonian constraint equation for ^' 
affects the diagonal elements of the three-metric produc- 
ing a nonlinear even-parity distortion. If both Hp and 
He vanish, undistorted Schwarzschild results. If they are 
present, they generate a linear odd-parity perturbation 
directly through the extrinsic curvature, and a second or- 
der even-parity perturbation through the conformal fac- 
tor Vt. Hence, the system will have mixed odd- and even- 
parity distortions mixed together at different perturba- 
tive orders. As we will see below, because even- and 
odd-parity components are cleanly separated in this way, 
and the background geometry is explicitly Schwarzschild, 
it is straightforward to construct analytic, linearized ini- 
tial data for these distorted black holes, which can then 
be evolved with the perturbation equations. 

In summary, our initial data sets contain both even- 
and odd-parity distortions of a Schwarzschild black hole, 
and are characterized by parameters (Qoj '^'j %j f )i where 
Qq determines the amplitude of the distortion,/!' deter- 
mines the angular pattern, 770 determines the radial lo- 
cation (with ryo = being the black hole throat), and a 
determines the radial extent of the distortion. For sim- 
plicity of discussion, all cases we will consider in this 
paper have the form (Qoj ?yo = 2, cr = 1). 



B. Two Perturbation Formalisms 



1. Metric Perturbations 



of Regge- Wheeler and Zerilli, is given in terms of the 
three-geometry metric perturbations. We will use the 
Moncrief formalism here as already described in Paper I. 

For special combinations of the perturbation equa- 
tions, a wave equation, the famous Regge- Wheeler equa- 
tion, resulted for a single function (p^im)'- 
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9t2 9r*2 
Here r* = r + 2M ln(r/2M — 1), and the potential 
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Because we are considering only axisymmetric perturba- 
tions all components with m ^ vanish identically. We 
will subsequently suppress the m labels. 

Moncrief showed that one can define a gauge-invariant 
function, that is invariant under infinitesimal coordinate 
transformations (gauge transformations), which is de- 
fined for any gauge via 
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(11) 



which satisfies the Regge- Wheeler equation above. As 
the Regge- Wheeler equation is a wave equation, in order 
to evolve the function (f> we must also provide its first 
time derivative, which is computed then directly through 
the definition of the extrinsic curvature of the perturbed 
Schwarzschild background: 



f^A. 2 / 2M 
OfCpi = — I 1 
r 



1 _ ^t-Wodd 
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This general prescription of linear Schwarzschild per- 
turbations simplifies dramatically in the present case. As 
discussed in Sec. HA above, the three-metric contains 



only even-parity perturbations, and those appear only at 
second order. Hence, for linearized treatment both the 
even- and odd-parity perturbation functions vanish in the 
initial data! To first order the metric is described by 
the Schwarzschild background and perturbed initial data 
consists solely of odd-parity extrinsic curvature contri- 
butions. Even-parity modes appear only at higher order, 
and are not considered in our comparisons here. 

For the specific initial data given in the previous sec- 
tion we obtain 



The theory of metric perturbations around a 
Schwarzschild hole was originally derived by Regge and 
Wheeler |2^ for odd-parity perturbations and by Zer- 
illi [ p7[ for even-parity ones. The spherically symmetric 
background allows for a multipole decomposition even 
in the time domain. Moncrief has given a gauge- 
invariant formulation of the problem, which like the work 
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n'=3,/=3, Q„=1,M^,„=1,Ti„=2,o=l 
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FIG. 1. The initial data for the Moncrief variable. For 
n' — 3 the only linear content is the i — 3 multipole. The 
initial value of (j> vanishes for linear odd-parity perturba- 
tions, since our choice of initial data only allows for second 
order even-parity perturbations of the initial three-metric. 
Our choice of the initial extrinsic curvature generates an 
almost Gaussian dt(j} that sits near the maximum of the 
Regge- Wheeler potential. 



where for n' = 3, £ = 3 and n' = 5, = 3, 5 the numerical 
coefficients coming from the multipole decomposition of 
the extrinsic curvature are 
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Although these initial data could be obtained numerically 
via an extraction process described in Paper I ||l^, it is 
not necessary to do so in this case with a clear analytic 
linearization. 

In Fig. |l| we plot an example of these analytic initial 
data. We are now ready to evolve linearly these data 
with the Regge- Wheeler equation. 



2. Curvature Perturbations 
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where M is the mass of the black hole, a its angu- 
lar momentum per unit mass, E = -f cos^ 9, and 
A = — 2 Mr -I- a^. The source term T is built up 
from the energy- momentum tensor p9[ |. Gravitational 
perturbations s = ±2 are compactly described in terms 
of contractions of the Weyl tensor with a null tetrad, 
which components (also given in Ref. pof) conveniently 
chosen along the repeated principal null directions of the 
background spacetime (Kinnersley choice) 



ip{t,r,9,(p) = 



P 1pi = -p Cnmnm for S = 
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for s — 



(19) 



where an overbar means complex conjugation and p = 
\/{r — ia cos 9) . This field represents either the outgoing 
radiative part of the perturbed Weyl tensor, (s = —2), 
or the ingoing radiative part, (s = +2). 

For the applications in this paper we will consider s — 
—2, since we will study emitted gravitational radiation 
and a = 0, i.e. perturbations around a Schwarzschild 
black hole. In general, for the rotating case, it is not 
possible to make a multipole decomposition of ■0 which 
is preserved in time. So to keep generality we shall use 
Eqs. (3.1) and (3.2) of Ref. |§ to build up our initial ^^4 
and (3tV'4 not decompose into i multipoles. The analytic 
expressions for the distorted black hole initial data sets 
considered in this paper, are: 
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sin^ 9[A^ 1 - 2^dr,qG - 29^'<Zg]} (20) 



There is an independent formulation of the perturba- 
tion problem derived from the Newman-Penrose formal- 
ism py] that is valid for perturbations of rotating black 
holes. This formulation fully exploits the null structure 
of black holes to decouple the perturbation equations into 
a single wave equation that, in Boyer-Lindquist coordi- 
nates {t, r, 9, ip), can be written as: 
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n'=3, Q„=1,M^„„=1,ii„=2,a=l, 9=71/4 
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FIG. 2. The initial data for the Weyl scalar ^/)4 and its 
time derivative as given in Eq. (po|). The general form of ip4 
resembles that of the Gaussian-like dt4>t- Note that 9tV'4 is 
nonvanishing in contrast to the initial (l>e ~ 



+ (3 + 2n' - n'^)r]d,,qG + [{n'^ + n' - 2) sin^ 6 
-(n'2-2n'-3)](5M-2r)r gc}- (21) 

As expected for pure odd-parity perturbations, only 
the imaginary part of "04 is nonvanishing. Also note that 
unlike (f>( — 0, 9t^4 does not vanish initially. This is 
because ip4 and its time derivative depend on both, the 
3-geometry and the extrinsic curvature while the Mon- 
crief function only depends on the perturbed 3-geometry 
and its time derivative only on the extrinsic curvature. 
We plot the initial data in Fig. ^. We evolve of these ini- 
tial data via Teukolsky cquation(^) using the numerical 
method is described in Ref. 

Since, after all, we are computing perturbations on a 
Schwarzschild (a = 0) background there must be a way 
to relate the metric and curvature approaches. In fact 
the relations between ■04 and Moncrief even- and odd- 
parity waveforms in the time domain have been found in 
Ref. and tested for the even-parity case in Ref. [p3| * 

Here we can perform the same kind of cross check for 
the odd parity modes. From the equations in section II. B 
of Ref. Q or (2.9) in Ref. Q we obtain a relation that 
holds at all times 



'{i + 2y. 

(^-2)! 



*Note that this relation among waveforms is only valid at 
first perturbative order. When nonlinearities are included the 
two approaches may give widely different results [p3[. 
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and that we can integrate to give us '04 from 0^ evolved 
with the Regge- Wheeler (^ equation instead of evolving 
V'4 directly with the Teukolsky equation (|l8|). 



C. Axisymmetric Nonlinear Evolutions 

The 2D fully nonlinear evolutions have been performed 
with a code, Magor, designed to evolve axisymmetric, 
rotating, highly distorted black holes, as described in 
Ref. ||8|-[l0|]. Magor has also been modified to include 
matter flows accreting onto black holes |3^ , but here we 
consider only the vacuum case. 

In a nutshell, this nonlinear code solves the complete 
set of Einstein equations, in axisymmetry, with maximal 
slicing, for a rotating black hole. The code is written in a 
spherical-polar coordinate system, with the rescaled ra- 
dial coordinate rj that vanishes on the black hole throat. 
An isometry operator is used to provide boundary con- 
ditions on the throat of the black hole. All three com- 
ponents of a shift vector are employed to keep all off 
diagonal components of the metric zero, except for the 
ge,^ component, which carries information about the odd- 
parity polarization of the radiation. For complete details 
of the nonlinear code, please see Ref s. &]IO| , p^ . 

The initial data described in Sec. II A above are pro- 



vided through a fully nonlinear, numerical solution to the 
Hamiltonian constraint. The code is able to evolve such 
data sets for time scales of roughly t < lO^M, and study 
such physics as horizons and gravitational wave emission. 

Consistently with the two different perturbations ap- 
proaches there arc two methods we use to extract infor- 
mation about the gravitational waves emitted during the 
fully nonlinear simulation: metric based gauge-invariant 
waveform extraction and direct evaluation of the curva- 
ture based Newman-Pcnrose quantities, such as ipi. 

The first method has been developed and refined over 
the years [ p5| , p6p^j37| , |38| , ^ to compute waveforms from 
the numerically evolved metric. Surface integrals of var- 
ious metric quantities are combined to build up the per- 
turbatively gauge invariant odd-parity Moncrief func- 
tions. These can then be compared directly with the 
perturbative results. 

A second method for wave extraction is provided by 
the calculation of the Weyl scalar -04 which is coordinate 
invariant but depends on a choice of tetrad basis. For our 
numerical extractions we follow the method proposed in 
Ref. To define the tetrad of their form numerically 
we thus align the real vector (which can be thought of as 
providing the spatial components of l^ and l^) with the 
radial direction. The complex vectors and point 
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within the spherical 2-surface. At each step, a Gram- 
Schmidt procedure is used to ensure that the triad re- 
mains ortho-normal. The tetrad assumed by this method 
is not directly consistent with the one assumed in the per- 
turbative calculation, but for the a = it can be made 
consistent by a type III (boost) null rotation which fixes 
the relative normalization of the two real-valued vectors. 
We have found that the transformation np —> A~^lo and 
Ip Alo where A — ^2/{l — 2M /r) fixes the normal- 
ization appropriately. For the general a ^ case this 
would be insufficient, and we would instead use the more 
general method proposed in Ref. | |3C| ]. 

D. Full 3D Evolutions with Cactus 

The last of our approaches for evolving these distorted 
black hole data sets utilizes full 3D nonlinear numerical 
relativity, and is based on Cactus. More general than 
a numerical code, the Cactus Computational Toolkit is 
actually a general parallel framework for numerical rela- 
tivity (and other sets of PDE's), that allows users from 
various simulation communities to gain high performance 
parallelism on many platforms, access a variety of compu- 
tational science tools, and to share modules of different 
evolution methods, initial data, analysis routines, etc. 
For the relativity community, an extensive suite of nu- 
merical relativity modules (or thorns in the language of 
Cactus) is available^, including black hole and other ini- 
tial data, slicing routines, horizon finders, radiation indi- 
cators, evolution modules, etc. 

For this paper. Cactus was used to assemble a set of 
3D initial data, evolution modules, and analysis routines 
needed for the comparisons with Magor and the two per- 
turbative approaches described above. All operations 
have been carried out in 3D Cartesian coordinates, from 
initial data to evolution to waveform extraction. The 
initial data are computed as in the Magor code, in a 
polar-spherical type coordinate system, and interpolated 
onto the Cartesian coordinate system as described in Pa- 
per I. The evolutions are carried out with a formulation 
of Einstein's equations based on the conformal, trace- 
free approach developed originally by Shibata and Naka- 
mura and Baumgarte and Shapiro and further 
tested and developed by Alcubierre, et. al., as described 
in . Due to certain symmetries in these initial data 

sets, the evolutions can be carried out in an octant, in 
Cartesian coordinates. However, we have chosen in this 
case to use the full 3D Cartesian grid, as enough mem- 
ory is now available to run sufficiently large scale evolu- 



Q„=2, n'=3,/=3, Ro,,=15M 



0.25 



0.15 



-0.15 



-0.25 





A 

/ \ / v/"---.^-^ 




MAGOR (300x59) 




1 Regge-Wheeler 



0.0 20.0 40.0 60.0 80.0 



m 

FIG. 3. The £ = 3, odd-parity Moncrief waveform, ex- 
tracted from the fully nonlinear 2D evolution code (solid line) 
and compared to the fully linear evolution of the same data, 
obtained by the Regge-Wheeler equation, for the low ampli- 
tude case {Qo ~ 2,n' = 3,770 = 2, a = 1). The excellent 
agreement shows that both approaches are accurate, and that 
the black hole physics is operating in the linear regime. 

tions that cover the entire spacetimc domain of interest, 
as would also be necessary when considering the general 
black hole inspiral problem. For comparison with the 2D 
code we extract waveforms via the same gauge-invariant 
Moncrief approach. In this case, surface integrals are 
carried out on the Cartesian based system by coordinate 
transformations and interpolation onto a coordinate 2- 
sphere, as described in Ref. |l9j . 

Further details of the individual simulation parameters 
are provided as needed when discussing the results below. 

III. RESULTS 

Here we compare the results of evolving the odd-parity 
distorted black holes by the four techniques described 
above. We consider two classes of distortions (n' = 3 and 

= 5) with different angular distributions, and various 
amplitudes to include cases of linear and distinctly non- 
linear dynamics. For the n' = 3 case the distortion is 
(linearly) pure £ = 3, while the n' — 5 case encodes a 
mix of £ = 3 and i — 5 distortions in the initial data. 



A. Comparison of Nonlinear Evolutions with 
Regge-Wheeler Theory 



t Cactus is well doc umented and can be downlo aded freely 



http: / /www. cactuscode.org . For more 



from a web server at 
information on Cactus, its use in numerical relativity and 
other fields, please see the web pages. 



In this subsection we compare the 2D nonlinear 
(Magor) evolutions with the results of the Regge-Wheeler 
perturbative approach. We first consider the nonlin- 
ear evolution of a family of data sets with parameters 
{Qo,n' = 3,770 = 2,0- = 1). 
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FIG. 4. The £ = 3, odd-parity Moncrief waveform, 
extracted from the fully nonlinear 2D evolution code, 
for a series of evolutions with initial data parameters 
{Qo,n = 3, 770 = 2, (J = 1). The waveforms are extracted at 
isotropic coordinates f = 15M. For n' = 3 we only have I = ?> 
linear contributions. We normalized waveforms by the ampli- 
tude Qo in order to study the linear and nonlinear regimes. It 
is observed that for Qo < 8 the linear regime is maintained, 
while for Qo = 32 nonlinearities are well noticeable. The ef- 
fect of nonlinear contributions increases the scaled amplitude 
of the waveform and increases its frequency. This indicates 
that the final ringing black hole is significantly smaller than 
the initial mass of the system. 

For low amplitude cases with Qo < 8, we are in the 
linear regime and even the nonlinear evolutions exhibit 
strongly linear dynamics. In Fig. ^ we show t — i wave- 
form results obtained from the 2D nonlinear code, for 
the case Qo = 2, and compare with Regge- Wheeler evo- 
lutions of the {(j),dt4') system. The agreement is so close 
that the curves cannot be distinguished in the plot. The 
perturbative-numerical agreement is equally good with 
the other linear waveforms at low amplitude so we will 
leave the perturbative results out of the plots and focus 
on the transition to nonlinear dynamics. 

In Fig. 1^ we show the I — i gauge-invariant Moncrief 
waveforms for a sequence of such evolutions of increasing 
amplitude Qo- The waveforms have all been normalized 
by the amplitude factor Qo — Qo /Af ^ to accentuate non- 
linear effects. If the system is in the linear regime, the 
normalized waveforms will all line up, as is clearly the 
case in the regime Qo < 8. For the large amplitude case 
Qo — 32, the normalized waveform is much larger, indi- 
cating that here we are well into the nonlinear regime. 

We now consider the transition to nonlinear dynam- 
ics in our second family of data sets, given by parame- 
ters (Qo,^^' = 5,7/0 = 2,cr = 1). These data sets have 
a linear admixture of both £ = 3 and € = 5 perturba- 
tions, and should contain waveforms of both types. In 
Fig. B, we show the results of the fully nonlinear evolu- 
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FIG. 5. We show the £ — 3 normalized odd-parity Mon- 
crief waveform for {Qo,n' — 5,rio — 2,a — 1) initial data, 
extracted from the fully nonlinear 2D code for a variety of 
amplitudes Qo. Again, the system is clearly linear for Qo < 8 
and nonlinearities cause an increase in amplitude and fre- 
quency of the wave. 

tion with the Magor code in maximal slicing, extracting 
the £ = 3 gauge-invariant Moncrief waveform for vari- 
ous amplitudes. The waveforms are again normalized by 
the amplitude factor Qo. In this case the nonlinearity is 
somewhat weaker at Qo ~ 32 so we have included the 
Qo = 64 curve in the figure. For the £ = 5 case, shown 
in Fig. ^, the higher frequency of the quasi-normal ring- 
ing makes it easier to appreciate the nonlinearities at 
Qo = 32. The plots indicate that again the dynamics are 
quite linear below Qo = 8 

The waveforms we have shown so far are the only ones 
predicted to linear order in perturbation theory. We 
would need to apply higher order perturbation theory 
to predict waveforms for the even-parity or higher-^ odd 
parity modes. Nevertheless general considerations from 
the perturbative point of view do provide some expec- 
tations on the scaling of the other nonlinear waveform 
modes within the families considered here. We return to 
the n' — 3 family for an example. The leading contri- 
bution to the n' = 3, £ = 5 odd-parity waveform comes 
from the cubic coupling of the first order £ = 3 odd- 
parity mode discussed above (including the coupling of 
the £ — 3 odd-parity mode with the second order even- 
parity £ — 2 mode expected via the source term contribu- 
tion to the solution of the Hamiltonian constraint in the 
initial data. Thus, this wave component should appear at 
the third perturbative order. We verify this expectation 
by plotting the numerical results for the £ = 5 odd-parity 
waveforms scaled this time by Qo 

Fig0. Although the 
magnitudes of these waveforms are far smaller than those 
of the £ = 3 mode we again see very nice agreement, be- 
low Qo — 8, with the perturbative expectation, that the 
waveforms should superpose. 
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FIG. 6. We show the £ — 5 normahzed odd-parity Mon- 
crief waveform for {QQ,n' = 5, ryo = 2, cr = 1) initial data, 
extracted from the fully nonlinear 2D code for a variety of 
amplitudes Qo. The regime is clearly linear for Qo = 2 and 
nonUnear components appear for Qo = 32. The increase in 
frequency and amplitude of the wave are seen here also. 



Let US now take another look at Figs, ^-jg to con- 
sider what is happening as we move into the nonlinear 
regime where the waveforms no longer superpose. In all 
the graphs we see the same general features arising as be- 
gin to drive the system into the nonlinear regime. These 
are higher frequency ringing, and larger amplitudes for 
the later parts of the waveform. At Qo = 32 the n' = 3 
case shows a roughly 10% increase in frequency compared 
to 5% for the n' = 5 case. Since the final state of the 
system will be a black hole, we expect quasi-normal ring- 
ing in the late-time behavior of the system regardless of 
the size of the initial perturbations. This is indeed what 
we see in the waveforms, except that the ringing is at a 
higher frequency (relative to the initial mass of the sys- 
tem) than we expected. This indicates that the mass of 
the final ringing black hole has less mass than the ADM 
mass of the initial data. The perturbations have grown 
large enough to generate radiation amounting to a no- 
ticeable fraction of the total ADM mass leaving behind 
a slightly smaller black hole. The smaller mass of the fi- 
nal black hole is also consistent with a larger amplitudes, 
since the scaled perturbation Qq = Qq/NP is larger rel- 
ative to a smaller mass black hole. The arrival time of 
the wave pulse is not strongly affected by the change in 
mass because the time and wave extraction points are 
both scaled against the initial mass. 



B. Comparison of Nonlinear Evolutions with 
Teukolsky Theory 
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FIG. 7. The Moncrief waveform for a purely nonlinear 
mode. For n' = 3, the £ = 5 multipole is generated by cubic 
products of the odd-parity wave (squares generate even-parity 
ones). Accordingly we normalized waveforms by Qq. Still 
higher nonlinearities switch on for Qo — 32 and show the 
generic increase in the frequency of the wave. 



We now turn to the curvature based Teukolsky ap- 
proach to perturbative evolution for black hole space- 
times. As motivated above, this is a much more power- 
ful approach that will enable perturbative evolutions of 
both polarizations of the gravitational wave, and evolu- 
tions of distorted black holes with angular momentum, 
and without the need for multipolar expansions. The key 
difference for analysis within this (Teukolsky) formalism 
is that , in the general case, we no longer have the ben- 
efit of a time-independent separation into multipolcs. In 
this first step towards the transition between the met- 
ric perturbation (Moncrief) approach and the curvature 
based (Teukolsky) approach, we consider the same data 
sets studied above, which have linear perturbations only 
for odd-parity, nonrotating black holes. We will consider 
more general systems in future papers. 

We now consider evolutions of the distorted black hole 
data set (Qo — = 3,770 = 2,cr = 1). The per- 

turbative initial data for -04 and dttpi, needed for use 
in the Teukols ky ev olution, have been obtained as de- 
scribed in Sec. II B 2| above. These data are then evolved 
and recorded at the same coordinate location as before 
(r = 15A/) for comparison with the previous results. In 
Fig. ^ we show the results of the full 2D nonlinear evolu- 
tion, obtained with Magor in maximal slicing, the Teukol- 
sky evolution, and the Regge- Wheeler evolution obtained 
previously. The solid line shows the results of V'4 obtained 
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FIG. 8. The Weyl scalar tp4 as seen by an observer located 
at isotropic coordinates f — 15M and — 7r/4. For n' = 3 we 
have a pure £ — 3 linear contribution although no multipole 
decomposition was made of tp4. We compare here the results 
of integrating full (axially symmetric) Einstein equations nu- 
merically with the program Magor using a resolution of 300 
radial and 59 angular zones with the linear evolution of initial 
data via the Teukolsky equation or the Regge- Wheeler equa- 
tion and then use the transformation equations ( |22| to build 
up ip4 [only possible in the nonrotating case]. 



with the Teukolsky equation, observed at a constant an- 
gular location 6 — 7r/4, and the dotted Hne shows the 
result of the Magor evolution at the same location, with 
■04 extracted from the full nonlinear simulation as de- 
scribed in Sec. II. C. 2 above. The results agree extremely 
well except at very late times, when the nonhnear results 
are affected slightly by coarse resolution in the outer re- 
gions of the numerical grid. We also verify here that 
the results of the Regge- Wheeler evolution, transformed 
to provide the same function ip4 according to Eq. ( |2^ ) 
above, agree with the Teukolsky evolution. The results 
of the two perturbative approaches are indistinguishable 
in the graph. 

We now examine the other family of distorted black 
holes with the choice of angular parameter n' = 5. The 
initial data were obtained as before, and evolved with 
the nonlinear Magor code and the Teukolsky code. The 
results are shown in Fig. ^, where we see excellent agree- 
ment between the two plots. But notice that the wave- 
form does not show the clear quasi-normal mode appear- 
ance that one is accustomed to in such plots. This is 
because this data set has a roughly equal admixture of 
both i = 3 and £ — 5 components of radiation, and the 
curvature based -04 approach is not decomposed into sep- 
arate multipoles. This waveform shows a clear beat of the 
two i — S and i — 5 components. 



FIG. 9. The Weyl scalar ■i/'4 for the n' — 5 initial data. We 
observe the beating of the £ — 3 and £ = 5 components since 
we are not making any multipole decomposition. 



C. Comparison of nonlinear 2D and full 3D codes 

Having successfully tested the 2D fully nonlinear code 
Magor for odd-parity distortions against perturbative 
evolutions, we can now test the 3D code Cactus against 
the 2D one . In Cactus, using the same procedures as 
in Sec. |[I C] , the initial data are evolved in the full (no 
octant) 3D mode, with a second order convergence al- 
gorithm, maximal slicing, and static boundary condi- 
tions. Note we perform the conformal-traceless scheme 
1^,^ for this evolution. The first observation is that 
we have to solve the initial value problem taking into ac- 
count all nonlinearities, even if we are in the linear regime 
{Qo < 10), since small violations of the Hamiltonian con- 
straint contaminate the outgoing waveforms. 

The runs presented in Figs. |l^, |l^, and ^ show 
very nice agreement with the 2D code (hence also with 
perturbation theory). Note that the spatial resolution 
(Ax^ = 0.15A/ — 0.3) is not high. Here we show wave- 
forms for t/M < 30. The runs do not crash afterwards, 
but become less accurate due to the low resolution and 
boundary effects, and even later to collapse of the lapse. 

The ^— modes shown in Figs. p^[|l^ are essentially dom- 
inated (for Qo = 2) by the linear initial distortion of the 
black hole. Those are the modes that we can compare 
with first order perturbation theory. Since we have two 
nonlinear codes we can now compare their predictions 
for modes dominated by nonlinear effects. That is the 
case of the odd mode £ — 5 when the initial data pa- 
rameter is n' = 3. This mode has a linear contribution 
only for ^ = 3. For ^ = 5 is easy to see that to gener- 
ate an odd mode we need at least cubic contributions. 
Thus this mode will scale as Qq. To be able to verify 
the agreement between the 2D and 3D codes we ampli- 
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FIG. 10. The £ — 3, odd-parity Moncrief waveform, ex- 
tracted from the fuUy nonhnear 3D evolution code Cactus 
(dotted line). The spatial grid consists of 256'^ points with 
a separation of 0.15M. The ADM mass of the black hole is 
M = 2.0 and (n' = 3, Qo = 2, t;o = 2, cr = 1). For compari- 
son we also plot the results of evolving the same initial data 
with the fully nonlinear 2D evolution code Magor (solid line) , 
which in turn has been tested against perturbation theory 
as shown in Fig. ^ Very good agreement is reached with a 
relatively low resolution of the 3D code. 
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FIG. 11. The 1 = 3, odd-parity Moncrief waveform pro- 
duced by the 3D code Cactus (dotted line). It corresponds to 
initial data with parameters (n' — 5,Qo = 2,?7o =2,a — 1). 
It also show very good agreement with the 2D results (solid 
line) which had been checked against perturbation theory as 
displayed in Fig. @. 



fied this mode taking Qq = 32 and checked the (almost) 
quadratic convergence of Cactus to the correct results as 
shown in Fig. |l^. 

IV. DISCUSSION 

We have completed a series of comparisons covering 
four different approaches for two classes of odd-parity 
distortions of Schwarzschild black holes. This includes 
2D and 3D nonlinear evolutions and for the first time, 
in both cases, a comparison of the odd-parity Regge- 
Wheeler-Moncrief formulation as well as the Teukolsky 
approach with numerical results. In all cases we find 
excellent agreement among the different approaches. We 
emphasize that these matchings have been achieved with- 
out the aid of any parameters and thereby stand as a 
strong verification of these techniques. 

Although the distorted black hole initial data config- 
urations we consider here are not necessarily astrophys- 
ically relevant, our analysis provides an example of the 
usefulness of perturbation theory as an interpretive tool 
for understanding the dynamics produced in fully nonlin- 
ear evolutions. In order to distinguish the cases of linear 
and nonlinear dynamics we simply show the output of the 
full nonlinear code, but we scale it by the factor Qq/M^ 
so that, if the system is responding linearly to Qq all the 
waveforms will lie exactly on top of one another. Using 
this procedure we are able to recognize the emergence of 
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FIG. 12. The £ = 5, odd-parity Moncrief waveform pro- 
duced by the 3D code Cactus (dotted line) with initial data 
having (n' = 5, Qo = 2, t^o = 2, cr = 1). Comparison with the 
2D results (solid line) show a very good agreement. See Fig. 
^ for the agreement between the 2D run and perturbation 
results. 



11 



n'=3, 1=5, Q„=32, R„,=9.22M, M,„=2.78 
0.06 I ' , ' , ' , — 




-0.06 I ' ' ' ' ' ' ' 1 

5.0 10.0 15.0 20.0 25.0 



t/M 

FIG. 13. The £ = 5, odd-parity Moncrief waveform pro- 
duced by the 3D code Cactus (dotted hne) with initial data 
having (n' = 3, Qo = 32,7;o — 2,a — 1) and AIadm ~ 2.777. 
This is a purely nonlinear mode, its leading term being cubic 
in the amplitude Qo- Comparison with the 2D results (solid 
line) show a good rate of convergence with the grid spacing 
(from Ax = 0.6M/2.777 to Ax = 0.3M/2.777. See Fig. Q for 
the purely 2D runs. 

nonlinear dynamics. Considering the mixing of pertur- 
bative modes also enables us to understand the results of 
one case which displays strictly nonlinear behavior, the 
£ = 5 waveform of the initial data with n' = 3 (see Figs. 
and |l|). This wave strictly vanishes to linear order 
in Qq and scales at lower amplitudes like Qq. The per- 
spective of perturbation theory allows us to create a full 
picture, identifying and explaining aspects of the nonlin- 
ear dynamics even when the perturbations are beyond 
the linear regime. In this case we find that linearized 
dynamics provide a very good approximation of the sys- 
tems' behavior until the radiation constitutes a signifi- 
cant portion of the initial mass, producing a smaller final 
black hole and, for example, higher quasi-normal ringing 
frequencies. 

Although we restrict to the case where the black hole 
system does not have net angular momentum, the ap- 
proach we develop in this paper is completely general, 
and can easily be extended to the case of distorted black 
holes with nonvanishing angular momentum. For this 
reason, we developed a procedure for using the Teukolsky 
equation to evolve the perturbations on a black hole back- 
ground, handling both the even- and odd-parity pertur- 
bations simultaneously, and providing the capability to 
deal with perturbations evolving on a Kerr background. 

In future work we expect to move in two directions: 
(a) We will apply the techniques developed here to the 
case of distorted, rotating black holes, to study nonhnear 
effects in the radiation of energy and angular momentum 
as well as to further develop the Teukolsky perturbative 
evolution paradigm for application to coalescing black 



hole initial data in a close-limit approximation, (b) We 
will use these techniques to evolve black hole systems, ei- 
ther from numerically generated initial data, or from par- 
tially evolved datasets that have reached a stage where 
they can be treated via perturbation theory. 
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